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ABSTRACT 

Using radiation magnetohydrodynamics simulations with realistic opacities and equa¬ 
tion of state, and zero net magnetic flux, we have explored thermodynamics in the 
inner part of protoplanetary discs where magnetic turbulence is expected. The ther¬ 
mal equilibrium curve consists of the upper, lower, and middle branches. The upper 
(lower) branch corresponds to hot (cool) and optically very (moderately) thick discs, 
respectively, while the middle branch is characterized by convective energy transport 
near the midplane. Convection is also the major energy transport process near the low 
surface density end of the upper branch. There, convective motion is fast with Mach 
numbers reaching > 0.01, and enhances both magnetic turbulence and cooling, rais¬ 
ing the ratio of vertically-integrated shear stress to vertically-integrated pressure by a 
factor of several. This convectively enhanced ratio seems a robust feature in accretion 
discs having an ionization transition. We have also examined causes of the S-shaped 
thermal equilibrium curve, as well as the thermal stability of the equilibrium solutions. 
Finally, we compared our results with the disc instability models used to explain FU 
Ori outbursts. Although the thermal equilibrium curve in our results also exhibits 
bistability, the surface density contrast across the bistability is an order of magni¬ 
tude smaller, and the stress-to-pressure ratios in both upper and lower branches are 
two orders of magnitude greater, than those favored in the disc instability models. It 
therefore appears likely that FU Ori outbursts are not due solely to a thermal-viscous 
limit cycle resulting from accretion driven by local magnetic turbulence. 

Key words: accretion, accretion discs — instabilities — magnetohydrodynamics 
(MHD) — protoplanetary discs — radiative transfer — turbulence 


1 INTRODUCTION 


Protoplanetary discs are accretion discs of mass accretion 
rates ranging from 2 x 10 -12 to 5 X lO _8 M 0 yr _1 and above 
Herczeg & Hillenbrand 2008). A most plausible mech¬ 


anism for driving accretion discs is magnetic turbulence 
caused by magneto-rotational instability (MRI) (a compre¬ 
hensive review by |Balbus fe Hawley|1998] and references are 
therein). However, the ionization fraction of gas is generally 
low for most of radii in protoplanetary discs since thermal 
ionization doesn’t work due to low temperatures. There¬ 
fore, non-ideal magnetohydrodynamics (MHD) effects are 
likely to suppress MRI there (Sano & Miyama 1999; Sano & 


|5h 


2002[ [20031 recent work s by, fo r example, |Bai Stone 


(2011); Wardle &; Salmeron (2011); Kunz &; Lesur (2013)) 


although non-thermal ionization by cosmic rays and the stel¬ 


lar X-rays may revive MRI in some limited places (Gammie 
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1996 |Sano et al.|2000; Ilgner & Nelson|2006||Turner &; Sano| 


2008). Pure hydrodynamic instabilities are also considered 


for such MRI-dead zones (Nelson, Gressel, & Umurhan 2013 
Klahr & Hubbard 2014). 


On the other hand, MRI must play an important role on 
driving the mass accretion in the inner region where thermal 
ionization is expected to work (Ar mitage| 2011). To evaluate 
the effect of thermal ionization on MRI turbulence, correct 
thermodynamics (Hirose & Turner 2011 Flock et al.|2013 ) is 
crucial since thermal ionization has strong temperature de¬ 
pendence. Thermodynamics has another significance in the 
inner region. It may regulate the final accretion onto the cen¬ 
tral star since the local relation between the surface density 
X and the mass accretion rate M is uniquely determined for 
the disc in thermal equilibrium. The relationship M = M(X) 
is called a thermal equilibrium curve, which is obtained by 
solving the vertical structure of the disc. Of particular inter¬ 
est is when the gas temperatures are around the hydrogen 
ionization temperature ~ 10 4 K, where a thermal equilib- 
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rium curve generally exhibits bistability. Then, disk annuli 
are expected to follow a limit-cycle, switching episodically 
between high and low accretion states. Actually, episodic 
high accretion states of M ~ 10 -4 Moyr -1 , called FU Ori 
outbursts, are observed in protoplanetary discs, and some 
authors applied the (hydrogen ionization) disc instability 
model (DIM) to explain them (|Kawazoe & Mineshige|1993;; 
Bell fe Lin|1994| ). 


The DIMs were originally developed to explain dwarf 
nova outbursts (Osaki|p~974| |HoshI |1979| Meyer &; Meyer- 


iofmeister|1981 Cannizzo et al.|1982[|Faulkner et aM1983| 


vlineshige & Osaki 1983, for a recent review, see Lasota 


2001)). They employ the a prescription of the shear stress 
Shakura &; Sunyaev|1973 ) and the mixing length theory for 


thermal convection. Consequently, in the DIMs, a thermal 
equilibrium curve strongly depends on two free parameters, 
a , the ratio of vertically-integrated shear stress to vertically- 
integrated thermal pressure, and the mixing length in terms 
of the disc scale height, both of which need to be tuned. On 


the other hand, Hirose et al. (2014) (hereafter Paper I) have 


derived a thermal equilibrium curve for dwarf nova discs 
from the first principles using radiation MHD simulations 
with realistic mean opacities and equation of state (EOS). 
They have also found that a is enhanced by a factor of 
several due to convection near the low surface density end of 
the upper branch of the thermal equilibrium curve. That is 
consistent with an observational fact that the high accretion 
state in dwarf novae characteristically exhibits a of order 0.1 
( Smak||1999 King et aL||2007 ). 


In this paper, using the same methods as in Paper I, 
we explore thermodynamics in the inner region of proto¬ 
planetary discs where MRI turbulence is expected to drive 
accretion. We also discuss whether FU Ori outbursts can be 
attributed to the thermal-viscous limit cycle in the frame¬ 
work of the DIMs. Since we focus on the intrinsic thermo¬ 
dynamics of accretion discs, we don’t consider the stellar 
irradiation here. Also, we don’t include, for simplicity, non¬ 
ideal MHD effects and the net vertical magnetic flux. These 
are planned to be treated in the future works. 

This paper is organized as follows. After we briefly de¬ 
scribe our methods in Section [2j we present the thermal 
equilibrium curve and the basic thermodynamical proper¬ 
ties in Section [ 3 ] In section^] we shed light on some points 
regarding the thermal stability that were not discussed in 
Paper I. The issue of FU Ori outbursts is also discussed in 
the section. Finally, we give summary in Section [5] 


2 METHODS 

Here, we describe our methods briefly since we employ the 
same methods as in Paper I and a complete description is 
given there. 


2.1 Basic Equations 


The basic equations are ideal MHD equations and frequency- 
integrated moment equations of radiative transfer: 


9(P V ) , V7 K *P 


(1) 


+V-(H=-Vp+ r (VxB)xB+^F, 
ut 47T c 

( 2 ) 

De 

Q- t + V • (ev) = -(V • v)p - (4j tB(T) - cE) n P p, (3) 
<9 E 

+ V • (Ev) = -Vv : P + (47 tB(T) - cE) K P p - V • F, 


— -Vx(vxB) = 0, 


( 4 ) 

( 5 ) 


where p is the gas density, e the gas internal energy, p the 
gas pressure, T the gas temperature, E the radiation energy 
density, P the radiation pressure tensor, F the radiation 
energy flux, v the velocity field vector, B the magnetic field 
vector (in CGS emu units), B(T) — (JbT a /tt the Planck 
function (<tb, the Stefan-Boltzmann constant), and c the 
speed of light. We use a flux-limited diffusion approximation 
in the radiative transfer, where F and P are related to E as 
F = -(cA(R)/^rp)VF and P = f (R)E. Here A (R) = (2 + 
R)/(6 ± 3 R ± R 2 ) is a flux limiter with R = \VE\/(krpE), 
and f (R) = (1/2)(1 — f(R)) I ± (1/2)(3 — f (R))nn is the 
Eddington tensor with f(R) = A (R) ± A (R) 2 R 2 and n = 
VE/\VE\ ( [Turner fe Stone| 2001). 

The EOS p = p(e, T), Rosseland-mean opacity 
ftR (p, T), and Planck-mean opacity Kp(p, T) are tabulated 
beforehand, which are referred in the simulations. Since we 
assume no explicit resistivity (i.e. assume ideal MHD) and 
no viscosity in the basic equations, the turbulent dissipation 
occurs only numerically. The numerically dissipated energies 
are captured in the form of internal energy in the gas, effec¬ 
tively resulting in an additional term q^ iss in the gas energy 
equation ^ (Hirose et al. 2006). In contrast, the radiation 
damping, which is another dissipation mechanism, is fully 
resolved in the simulations ( Blaes et al.|201 1, See Appendix 
[A] for details). 


2.2 Numerical Methods 


We use the shearing box approximation to model a lo¬ 
cal patch of an accretion disc as a co-rotating Cartesian 
frame (x,y,z) with a linearized Keplerian shear flow v = 
— (3/2 )Qxy, where the x, y, and z directions correspond to 
the radial, azimuthal, and vertical directions, respectively, 


and y is the unit vector in the y direction (Hawley et al 


1995). In this approximation, the right hand side of the equa¬ 
tion of motion © has extra terms representing the inertial 
forces, —2 Qz x v + 3Q 2 xx — Q 2 zz , where x and z are the 
unit vectors in the x and z direction, respectively. Shearing- 
periodic, periodic, and outflow boundary conditions are used 
for the x , y , and z boundaries of the box, respectively (Hi- 
et al . 2006 ) 


The basic equations are solved time-explicitly by ZEUS 
using the Method of Characteristics-Constrained Transport 
(MoC-CT) algorithm except for the radiation-gas energy 
exchange terms ±(47 tB — cE)npp and the radiative diffu- 
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sion term —V • F. These terms are coupled and solved 


time-implicitly using Newton-Raphson iteration (Tomida et 


al. 2013) and the multi-grid method with a Gauss-Seidel 


smoother. 


2.3 Initial Conditions and Parameters 

A shearing box is characterized by a single parameter, the 
angular velocity Q, which appears in the inertial force terms 
and the shearing periodic boundary condition. We set Q to 


dwarf nova disc case in Paper I. The value corresponds to 
a distance of Ro = 5.89 x 10 11 cm = 0.0394 AU = 8.46R© 
from the central star of 1M©. Although was chosen arbi¬ 
trarily, the basic features of the thermal equilibrium curve 
will not depend on it strongly. (See, for example, Figure 3 
Kawazoe & Mineshige (1993}, where they show how the 


thermal equilibrium curve changes with Q.) Therefore, we 
expect that the radius we chose reasonably represents the 
inner region of protoplanet ary discs. 

The initial conditions for gas and radiation are specified 
by two parameters, the surface density Eo and the effective 
temperature T e ff 0 |j assuming the vertical hydrostatic and 
radiative equilibriums (see Appendix in Paper I). The initial 
magnetic field is a twisted azimuthal flux tube with zero net 
vertical flux placed at the centre of the simulation box. The 
field strength (typically 30 in terms of plasma beta) is chosen 
so that the initial development of MRI is resolved. The initial 
velocity field is the linearized Keplerian shear flow, whose x 
and z components are perturbed by random noise of 0.5 % 
of the local sound velocity. 

The box size is (1,4, 8) in terms of the initial disc thick¬ 
ness ho, and the numbers of cells are (32,64,256), in the x, 
y, and z directions, respectively. These values are the same 
as those employed in Paper I. We will discuss the numerical 
convergence of our results in Section |T6] Different box sizes 
are also used in some limited runs. For details of parameters 
of runs, see Table [l] 


3 RESULTS 


In this section, we present results of our simulations us¬ 
ing the same diagnostics in Paper I. The diagnostics are 
based on horizontally-averaged vertical profiles, which were 
recorded every 0.01 orbits. The horizontally-averaged verti¬ 
cal profile of quantity /, for example, is computed as 


/ave (-25 t) 


If dxdy 


( 6 ) 


where the integrations are done over the full extent of the 
box in x and y. 


3.1 Numerical Procedure 

Our procedure for obtaining a thermal equilibrium curve is 
as follows — (1) We start a simulation from the initial con¬ 
dition determined by a set of the two parameters (E 0 , T e ff 0 ). 
(2) We continue the simulation until the disc patch reaches 


1 Here, “0” denotes an initial value. 
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a thermal equilibrium that lasts at least for 100 orbits or 
experiences a thermal runaway. (3) For the former, we com¬ 
pute time-averaged effective temperature T e s and surface 
density E. Since a small amount of mass can leave a box 
from the vertical boundaries, E differs from Eo, typically by 
several percent (See Table [l}. — Repeating the above three 
steps for many sets of (Eo,T e ffo)> we obtain a thermal equi¬ 
librium curve, or the effective temperature as a function of 
the surface density, f e ff = T e ff(E). 

The time-averaged effective temperature and surface 
density in the above are computed as 


2<r B T e 4 ff = J (q ) dz, 

(7) 

E= f (p) dz. 

(8) 


The cooling rate per unit volume q is defined as 

-_dF z d((e + E)v z ) 

q = ~dT + - dz -’ (9) 

where F z is the z component of radiation energy flux F. 
Here and hereafter, the vertical integration f dz is done for 
the full extent of z. The brackets () denote time-averaging 
of a horizontally-averaged quantity; for example, 

(/> (z) = /ave(z, t)dt, (10) 

where averaging is done for a selected period of A in which 
the disc patch is in a quasi-steady state and the MRI near 
the midplane is reasonably resolved. The criterion (= 100 
orbits) for judging a thermal equilibrium of the disc patch 
was chosen arbitrarily, but, as shown in Table [l] the aver¬ 
aging period A is long enough (at least, several) in terms of 
the thermal time ttherm computed (a posteriori) as 


_ f (e + E) dz 
Uherm = f{q~)dz 


( 11 ) 


3.2 Thermal Equilibrium Curve 


Figure ^ shows the thermal equilibrium curve of the disc 
patch, T e ff = T e ff(E), consisting of the thermal equilibrium 
solutions defined in section |3.1| The right axis shows the 
corresponding mass accretion rate, 


Ff _ 4tt 2 gr B r e 4 ff 
“3 Q 2 ' 


( 12 ) 


The color indicates the time-averaged total Rosseland-mean 
optical thickness rtotai = f (pK r) dz. Table |T] lists time- 
averaged quantities for the runs shown in the figure. 

There are two distinct solution branches; the upper 
branch (? e ff > 5000 K, or M > 10 _5 M© yr -1 and E > 
6000 gem -2 = E m in) and the lower branch (T e ff < 1900 K, 
or M < 10 -7 M© yr -1 and E < 10000 gem -2 = E max ). For 
a limited range of surface density (E m m < E < E max ), there 
exist, for a single value of E, equilibrium solutions both on 
the two branches, exhibiting bistability. The bistability is 
caused by a strong temperature dependence of Rosseland- 
mean opacity around the hydrogen ionization temperature 
T ~ 10 4 K (e.g. Cannizzo 1993). That is, given a surface den¬ 


sity within the range, the disc patch can be hot and opaque 
(t total ~ 10 5 ' 5 ) on the upper branch or cool and less opaque 
(r to tai ~ 10 2 ' 5 ) on the lower branch. Here, the Kramers-type 
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opacities are responsible for the former while atomic and 


molecular opacities are for the latter (see also Figure 10 for 
the temperature dependence of opacities). Note that linearly 
unstable equilibrium solutions that appear in the DIMs (e.g. 
Kawazoe & Mineshige|l993) are excluded in our simulations 
where nonlinear evolution of the disc patch is calculated. 

The above features of the S-shaped thermal equilibrium 
curve are qualitatively similar to, but quantitatively differ¬ 
ent from the dwarf nova disc case in Paper I, in which 
is about 250 times larger than here. First, both the mini¬ 
mum effective temperature on the upper branch (~ 5000K) 
and the maximum effective temperature on the lower branch 
(~ 1900K) are about 30% smaller than those in Paper I. 
Such trend of the maximum/minimum effective tempera¬ 
tures with Q is also observed in the DIMs (e.g. Figure 3 in 
Kawazoe &; Mineshige|1993 ). Second, the bends of the ther¬ 
mal equilibrium curve occur around a surface density that is 
about 60 times larger (E ~ 10000 gem -2 ). This is required 
from the thermal equilibrium condition of a disc patch, 

2cr B T e 4 ff = J w V (f)dz = J pdz « ^QaT mid E, (13) 


where w r( j) is the shear stress and the a prescription is used 
in the second equality ( Cannizzo||19 93). From this relation 
(cr B T 4 ff « aT m id HE), E must increase by a factor of 250 x 
0.7 4 ~ 60 to compensate the increase of a factor of 250 in 
Q and the decrease of a factor of 0.7 in T e ff, provided that 
a and T m id don’t change much with Q, which is the case 
here. Third, E max /E m i n , or the relative range of E for the 
bistability, is smaller; it is ~ 3 in Paper I while it is ~ 1.6 
here. 

We see another solution branch, the middle branch 
(T e ff ~ 3000 K, E ~ 10000 gem -2 ), where the total opti¬ 
cal thickness (ftotai ~ 10 6 ) is even higher than on the upper 
branch. As shown in the following sections, what makes the 
solutions on this branch different from those having the same 
E on the upper branch is convection being the heat transport 
mechanism near the midplane. We observed thermal equi¬ 
libria for three different initial surface densities, Eq = 7949 


(ws0858), 8754 (ws0831 and wa0831), and 10542 (wb0828 
and wc0828) although the disc patch in ws0858 and wa0831 
collapsed eventually (indicated by downward triangles). Two 
runs, having the same Eo, slightly differ in the box size and 
the resolution (see Table [TJ . We will discuss the thermal 
stability of these solutions later in section |4~2] together with 
the two solutions at the high E end of the lower branch 
(ws850 and ws854), in which the disc patch flared in the 
end (indicated by upward triangles). 


3.3 Heat Transport 

Next, we examine heat transport mechanisms in the thermal 
equilibrium solutions. In Figure [2] time and horizontally- 
averaged vertical profiles of radiative heat flux T r - d (z), 
advective heat flux F adv(*)> and cumulative heating rate 
F^at (z) are shown for a typical upper-branch solution 
(ws0805), the solution at the low E end of the upper branch 
(ws0837), a middle-branch solution (ws0831), and a typical 
lower-branch solution (ws0800). Here, the heat fluxes are 


computed as 

F~M) = (F z ), (14) 

F- dv ^) = ((e + E)v z ), (15) 

FLtW = f «4ss> + (- (V • v)p) + <-V« : P» dz. 

Jo 

(16) 

From the energy equations © and Q, thermal energy 
balance in a steady state is written as 

9dlss -(V-u)p — Vu:P = V- F + V-((e + E)v) , (17) 

where the left hand side is the total heating rate (turbulent 
dissipation and compressional heating) and the right hand 
side is the total cooling rate (radiative diffusion and advec- 
tion). Note that the numerical dissipation rate q^ iss is not 
explicitly written in equation ©• Therefore, it is expected 
that 


Theat (^0 — T rad (z) + T adv (^), (18) 


which is the time and horizontally-averaged version of equa¬ 


tion Actually, equation (18) holds well in all solutions 
shown in the figure, confirming thermal equilibrium there. 
We note that the compressional heating — (V • v)p — X7v : P 
can act as virtual dissipation when radiative diffusion ex¬ 
ists. That is, radiation diffusion leads to photon damping 
and makes compressive heating irreversible. That is really 
observed in some solutions though it is a small fraction (10% 
or so) of the total dissipation when vertically-integrated (see 
Appendix [A|. The same mechanism, the radiation damping 
(Silk damping; Silk 1968), acts more efficiently in radiation- 
dominated discs ( Agol fe Krolik|1998 Blaes et al.||2011 ). 

Figure [2] also clearly shows that the main heat trans¬ 
port mechanism differs between the solutions. In the typ¬ 
ical upper-branch and lower-branch solutions (ws0805 and 
ws0800), radiative diffusion accounts for heat transport in 
the entire heights. However, in the solution at the low E 
end of the upper branch (ws0837) and the middle-branch 
solution (ws0831), it is advection that transports heat near 
the midplane. The heat advection is associated with thermal 
convection that is induced by high Rosseland-mean opaci¬ 
ties around T ~ 10 4 K. Actually, squared (hydrodynamic) 
Brunt-Vaisala frequency N 2 computed as 


N 2 


1 din (p) 


(ri) d\nz d\nz 


dln (p) _i n ( p ) dln ( r i> 


dlnz 


(19) 


is negative where advection transports heat. (Here, F± is the 
first generalized adiabatic exponent.) This, in turn, means 
that convection does not cancel the super-adiabaticity (Pa¬ 
per I). 

In Figure [3] snapshots on the x-z plane well contrast a 
“radiative” solution (ws0805, t — 71.6 orbits) and a “con¬ 
vective” solution (ws0837, t — 108.3 orbits). (These two so¬ 
lutions correspond to the top two panels in Figure [2] re¬ 
spectively.) In the convective solution, convective plumes, 
whose sizes are a fraction of the time-averaged pressure scale 
heigh© are clearly seen in the panels of specific entropy, 


2 Both the x and z axes in the figure are normalized with the 
time-averaged pressure scale height. 
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density and gas temperature. Also, upward (downward) con¬ 
vective plumes well correspond to upward (downward) ad- 
vective heat flux (ev z ), respectively. Note that the downward 
and upward convective plumes are highly asymmetric under 
the vertical stratification. The convective motion is strong 
enough to greatly deform the disc patch. In the radiative so¬ 
lution, on the other hand, the disc patch basically keeps its 
shape, only fluctuating because of the MRI turbulence and 
the spiral acoustic waves (Heinemann &; Papaloizou|[2009 ). 
There is no convection, as indicated by the entropy gradi¬ 
ent, and heat transport is done by radiative diffusion in the 
direction of the gradient of gas temperature. Note that the 
energy exchange between gas and radiation is so rapid that 
gas and radiation share a single temperature here. 


3.4 Vertical Profiles 


Figure^] shows time and horizontally-averaged vertical pro¬ 
files of density and pressures. Magnetic pressure is always 
dominant at high altitudes, supporting an extended expo¬ 
nential atmosphere (e.g. Hirose k, Turner||2011 ), while it is 
dominated by thermal pressure near the midplane. These 
are characteristic features generally seen in stratified MRI- 
turbulent discs. Note, however, that plasma beta near the 
midplane is rather small in the convective solutions; espe¬ 
cially, it is ~ 10 at the low E end of the upper branch 
(ws0837), which is about ten times smaller than the char¬ 
acteristic value. Also, in the convective solutions, density 
is peculiarly flat near the midplane. Such flat density pro¬ 
files are also observed in convective solutions in IBodo etl 


al. (2012), where they solved an energy equation with fi¬ 


nite thermal diffusivity and a perfect gas EOS. Radiation 
pressure is always smaller than gas and magnetic pressures 
except in the upper branch solution ws0805, where it is a 
fraction of gas pressure and larger than magnetic pressure 
near the midplane. 

Figure [5] shows time and horizontally-averaged vertical 
profiles of gas temperature T and the ionization fraction 
/ion- In the upper-branch solution ws0805, T > 10 4 K and 
/ion = 1 at almost entire heights. At the low E end of the 
upper branch (ws0837), such complete ionization still holds 
near the midplane, but T and f 1Q n drop to 3000 K and 10 -4 , 
respectively, beyond the photospheres. In the middle branch 
solution ws0831, T is just below 10 4 K and f lou is around 
10“ 1 near the midplane, and they drop to 2000 K and 10 —5 , 
respectively, beyond the photospheres. In the lower-branch 
solution ws0800, T is around 1000 K and roughly constant 
due to rather small optical thickness (as indicated by colors 
in Figure [TJ , and /i on is between 10“ 10 and 10 -8 . 


3.5 Enhancement of a 

Paper I found that the time-averaged a, computed as 

_ = Wxy = + {pv x (v v + 3Qx/2 )))dz 

“ P\ therm “ f(( P ) + (E) /3 )dz ’ ^ ' 

is enhanced near the low E end of the upper branch due to 
convection. Such enhancement of a is also observed in our 
simulations, suggesting that it doesn’t depend on £2 of the 
disc patch. 

Figure [6] (a) shows a as a function of E. For most of the 
© 0000 RAS, MNRAS 000, 000-000 


solutions, a is ~ 0.03, a typical value of the MRI turbulence. 
However, a increases as E decreases near the low E end 
of the upper branch, and exhibits a maximum value of ~ 
0.14 at the very end (ws0837). Also, middle-branch solutions 
show rather high values of a ~ 0.06. The direct reason for 
the high a on the upper branch can be seen in the trends 
of W xy and Ptherm with E (Figure [6] b). The stress W xy and 
pressure Ptherm correlate well in the solutions that show a 
of a typical MRI value. However, near the low E end of the 
upper branch, they deviate from the standard trend on the 
upper branch in opposite ways; as E decreases W xy goes 
up while Ptherm goes down, which makes a = W xy /Ptherm 
larger. Physically, the opposite trends of W xy and Ptherm 
can be attributed to convection as Paper I argued. First, 
W X y goes up because vertical convective plumes (Figure [ 3 ]) 
can create coherent vertical fields that seed the axisymmetric 
MRI. Second, Ptherm goes down because convection can cool 
the disc patch more than the increased stress, or dissipation 
heats it. In our case, the second effect is more prominent in 
enhancing a , judging from Figure [6] (b). 

Next, we examine the condition for convection to en¬ 
hance a in the MRI turbulence. Paper I introduced, as 
a barometer of convection, the advective fraction of heat 
transport computed as / a dv = ( /adv ), where 


/adv(t) 


J( F ady} Sgn(d{p t herm}<k: 

/ ((^adv) + (^rld)) Sgn(z){f>therm}dz' 


( 21 ) 


Here, ptherm = p + E/3 is used as a weight function, and the 
brackets {} denote a boxcar smoothing of a horizontally- 
averaged quantity for a single orbit. As shown in Figure [6] 
(c), solutions that show high a also show high / a dv, implying 
strong convection. However, it is not always true the other 
way around. For example, the values of / a d v are similar both 
near the low E end of the upper branch and on the middle 
branch, but a is smaller in the latter. Also, the high-E-end 
solutions on the lower branch show high / a d v while they show 
a of a typical MRI value. The discrepancy between a and 
/ a dv in the second example was also discussed in Paper I and 
they claimed that Mach number of the convective motion 
needs to be large for a to be enhanced. Here, we check their 
claim quantitatively by computing the Mach number of the 
convective motion as follows: 


1 f {(e +E)v z }sgn(z) 

f {ptherm}dz J {(e + E)}{c s } 


where the sound speed c s is defined as <f = 
(Fijp + 4/3-E7/3) //?. As shown in Figure [6] (d), M a d v actu- 
ally correlates with a much better than / a d v - Therefore, we 
confirm that convection alone is not responsible for high a , 
but its motion needs to be fast in terms of Mach numbers 
reaching > 0.01. 

Paper I reported that they observed convec¬ 
tive/radiative limit cycle near the low E end of the 
upper branch, where f&dv(t) switches between 0 and 1 
episodically. However, we didn’t observe such a clear 
“on/off” limit cycle in our simulations. For example, at the 
low E end of the upper branch (ws0837), convection occurs 
rather continuously with / a d v ~ 0.8 and the standard 
deviation of / a dv(t) is ~ 0.1. 

Although it is true that convection affects the dy- 
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namics of turbulence, MRI seems still primarily responsi¬ 
ble for driving the turbulence since the properties of the 
standard MRI turbulence are retained. As shown in Fig¬ 
ure [ 7 I W xy well correlates f (B 2 ) dz, and f (—B x B y ) dz = 
3 f (pv x (v y + 3Qx/2)) dz in all runs ( Sano et al.|2004 ). 


(2007) will be needed for them. However, such resolution 


3.6 Numerical Convergence 

Dependence on the grid resolution needs to be looked at 
when specific numerical values of a are quoted. For exam¬ 


ple, Fromang and Papaloizou (2007) have shown that the 


saturation of MRI turbulence in unstratified shearing box 
simulations without net vertical magnetic field and explicit 
dissipation depends on the grid resolution; a decreases as 
the resolution is increased. On the other hand, as for strat¬ 
ified shearing box simulations (of MRI turbulence without 


studies are highly computationally demanding for our cur¬ 
rent radiation MHD simulations, and thus are left for future 
works. We note that Paper I has done a limited resolution 
study for similar convective solutions to show that the en¬ 
hancement of a is not sensitive to the grid resolution. 


4 DISCUSSION 

4.1 Causes of the S-shape 

In the DIMs, a thermal equilibrium curve is generally S- 
shaped on the E—T e ff plane. That is, two stable solution 
branches of a positive slope are connected by an unstable 
solution branch of a negative slope. One way to see this is 


net vertical field and explicit dissipation) like ours, Hawley, to see equation (13), creT eff ~ oT m idDE, which says that 


Guan, & Krolik (2011) have proposed quantitative diagnos¬ 


tics by which the numerical convergence on the grid res¬ 
olution can be assessed. Specifically, the convergence met¬ 
rics that show clear resolution dependence and thus are ap¬ 
propriate for the assessment are the number of cells within 
the fastest-growing MRI wavelength in the z direction, the 
same but in the y direction, the ratio of magnetic stress 
—B x B y /A7 r to magnetic energy B 2 /Sn and the ratio of B 2 
to B 2 . Here, we compute the four convergence metrics as 

n - I VmW¥) (Ptherm) dz 2tT 

f( P t h erm)dz QAz’ (M) 

_ I \j (-®y) /** (p) (Ptherm) dz 2^ 

QV ~ f (Ptherm) dz nA^’ ^ 

_ f (( — BxBy/ 47r) / (B /8tt}) (ptherm) dz 
Qmag = fip^erm) dz 5 ^ j 

n _ / {(Bx) / (By)) (Ptherm) dz 

J mag — r , \ j 5 V ZD / 

J vPtherm) dz 

where the time and horizontally-averaged thermal pressure 
(ptherm) is used as a weight function^] 

The results are shown in Figure [8] According to |Haw-| 


ley, Guan, & Krolik (2011), numerical convergence is seen 
in those simulations that show Q z > 10 and Q y > 20, and 
ttmag — 0.3-0.4 and / mag > 0.15 are signatures of well- 
developed MRI turbulence. In our case, Q z and Q y meet 
the criteria of numerical convergence in most of the solu¬ 
tions. As for the signatures of well-developed MRI turbu¬ 
lence, c^mag meets the criterion while /mag is mostly ~ 0.11, 
about 25% smaller than the criterion. Looking at Figure 5 
in Hawley, Guan, &; Krolik (2011), however, 


u 


the resolution is increased, but begins to level off when it 
exceeds ~ 0.1. Therefore, we may conclude that MRI turbu¬ 
lence is marginally resolved in most of our solutions while it 
is perhaps underresolved in the four solutions on the lower 
branch that show / mag <0.1 and Q z ,Qy < 10. The values 
of a(~ 0.02) in the four solutions should be underestimated. 

The above convergence metrics may not be directly ap¬ 
plied to the solutions in which convection comes in, and di¬ 
rect resolution studies as done in |Fromang and Papaloizou] 


3 The results are almost unchanged when the time and 
horizontally-averaged density (p) is used as a weight function. 


the sign of the slope is determined by the dependence of 
T m id on T e ff. That is, E usually correlates positively with 
T e ff, but the correlation becomes negat ive when T m jd de¬ 
pends on T e ff more strongly than T 4 ff ( Cannizzof1993 ). The 
dependence of T m id on T e ff can be actually examined in our 
simulations. As shown in Figure [9] it is always weaker than 
T 4 ff on the solution branches, including the middle branch, 
while it becomes almost critical near the low E end of the 
upper branch. Therefore, our results also imply that the disc 
patch is thermally unstable when the dependence of T m id on 
T e ff is stronger than T e 4 ff . 

The dependence of T m id on T e ff can also be qualita¬ 
tively discussed with time and horizontally-averaged vertical 
profiles of the solutions ((T) (z), ( p) (z)) plotted on the T-p 
plane (Pojmanski 1986| ). As shown in Figure [To] two gaps di¬ 
vide a bunch of profile curves into three groups, which, from 
the left to the right, correspond to the lower branch, the 
middle branch, and the upper branch. The slope of a profile 
curve is always positive for all solutions, indicating that (T) 
and (p) are monotonically increasing from the surface to¬ 
ward the midplane. We don’t see any density inversion that 


is seen in some convective solutions in the DIMs (Pojman¬ 
ski 1986). The profile curves of the lower branch solutions 


are nearly vertical since the total optical thickness r to tai is 
relatively small. On the other hand, the profile curves of the 
solutions near the low E end of the upper branch are almost 
horizontal (p ~ const.) near the midplane due to convec¬ 
tion (as we saw in Figure [4| while they are almost verti¬ 
cal (T ~ const.) at high altitudes because of low opacities. 
Consequently, the gap between the upper-branch and the 
lower-branch profile curves is larger in T m id than in T e ff. In 
other words, a small change in T e ff results in a large change 
in T m id, or T m id strongly depends on T e ff between the up¬ 
per and lower branches. Note that the profile curves of the 
middle branch solutions are nearly parallel, indicating that 
T m id doesn’t strongly depend on T e ff (as seen in Figure [9]). 

Figure [To] also tells us about what is happening at the 
low (high) E end of the upper (lower) branch, respectively. 
Near the low E end of the upper branch, as E decreases, a 
convective region appears away from the midplane. The con¬ 
vective region finally extends to the midplane at the end of 
the branch, where the disc patch loses a thermal equilibrium 
( Pojmanski|1986' Cannizzo|1993| ). Note that the convective 
region on the upper branch well correlates with the region of 
high Rosseland-mean opacity (panel (b)), implying that the 
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convection is induced by the high opacities. On the other 
hand, near the high E end of the lower branch, convection 
appears at the midplane as E increases. Since the opacity is 
still low there (panel (b)), the main cause of convection may 
be low ri (panel (c)), which reduces the adiabatic temper¬ 
ature gradient. The disc patch seems to lose an equilibrium 
when the opacity source switches from atoms/molecules to 
H - at the end of the branch. 


4.2 Thermal Stability of Equilibrium Solutions 


When the cooling rate Q depends on temperature T more 
strongly than the heating rate Q + does, 


d In Q + ^ din Q 
dlnT < dlnT ’ 


(27) 


the disc patch is thermally (linearly) stable. The condition 
can be written as 


dlogQ 
d log Q+ 


(28) 


We examine this condition for selected thermal equilibrium 
solutions including those exhibit a thermal runaway in the 
end. Figure m shows time trajectories of the selected so¬ 
lutions on the logQ + -logQ _ plane. Here, Q + and Q~ are 
short-term average, vertically-integrated heating and cool¬ 
ing rates, 


Q+(t) = I {«,+„ - (V • v)p - Vi> : P} dz, (29) 

«'c>'/{f< 3 °> 


The centres of these trajectories are roughly on the line 
logQ - = logQ + , showing that the solutions are thermally- 
stable in a statistical sense. 

First, we see typical solutions on the upper branch 
(ws0803, the top panel) and on the lower branch (ws0800, 
the bottom panel), both of which are radiative solutions. If 
we look at the trajectory curve closely, it is almost horizon¬ 
tal on short timescales due to a rapid stochastic variation of 
Q + (£). Even in a statistical sense, the trajectory is roughly 
aligned with dlog Q - /dlog Q + = 1, rather than satisfying 
the condition (28). Therefore, the simple linear theory does 
not apply here. Next, we see the trajectories of the convec¬ 
tive solutions on the upper branch (ws0837) and the mid¬ 
dle branch (wa0831 and ws0831). They are more compact 
than those of the radiative solutions, and thus the slope is 
less clear. Also, it is not clear from their similar trajectories 
why the disc patch finally collapsed only in wa0831, and 
not in the other two cases. On the other hand, the trajec¬ 
tory of the solution at the high E end of the lower branch 
(ws0850), which also experienced a runaway finally, looks 
different from those mentioned in the above, and its slope is 
apparently less than unity, indicating thermal instability. 

The slope of the trajectory may be quantitatively eval¬ 
uated as JlogQ - /(51ogQ + , the ratio of the standard devia¬ 


tion of logQ (t) to that of logQ + (t). In Figure 12 we plot 


the ratio as a function of E to see the trend among the equi¬ 
librium solutions. The ratio is almost unity for most of the 
solutions (as we see for ws0803 and ws0800 in the above). 
However, it is obviously smaller for solutions near the ends 
of the upper and lower branches and on the middle branch. 


It is especially small (less than 0.7) for the solutions that 
finally experienced a thermal runaway (as indicated by tri¬ 
angles), which suggests that the ratio can be an indicator 
of “fragility” of the thermal e quilibrium (Hir ose et al.|2 009 ; 
Jiang, Stone, & Davis 2013). Other solutions with Q~ (t) 
varying less than the Q + (£) stayed in equilibrium up to the 
end of the simulations, but nevertheless are potentially un¬ 
stable since Q + (£) varies stochastically and Q~ (t) does not 
always respond on the same timescale. So, if Q + (t) wanders 
around the equilibrium or returns quickly after departing, 
then Q~(t) can respond and the system is stable. However, 
a sufficiently large and long-lasting excursion in Q + (£) builds 
up a heat excess or deficit so big that the patch of disk does 
not return to the equilibrium. Thanks to this fragility, the 
ends of the upper and lower branches are not well-defined 


(Latter & Papaloizou 2012). Another consequence is that 


the middle branch is long-lasting in some runs, but not in 
others with similar parameters, as also observed in Paper I. 

The stability of the middle branch solutions can be dis¬ 
cussed from another viewpoint. In the DIMs, it is explained 
that the portion of a negative slope of the S-shaped ther¬ 
mal equilibrium curve is thermally unstable du e to a strong 
temperature dependence of H - opacity (e.g. 


Kato, Fukue. 


&; Mineshige|2 008). That is, using ~ p 1 / 2 T g for H , the 


cooling rate anti-correlates with the midplane temperature 


Q~ 


rri4 

rp4 mid 

' 1 eff ~ 


—19/4 


T 

1 mid 


KrS E3/2 


(31) 


while Q + ~ ET m id assuming the a prescription, violating 
the thermally stable condition \27\ . This is, however, true 
when the disc patch is cooled only by radiative diffusion in 
which cTRTg ff « 4acT^ id /3^RE (a is the radiation constant). 
To put it the other way round, another cooling mechanism, 


or convection, may destroy the condition (31), stabilizing the 
disc patch. Actually, all the middle branch solutions have ri 
of the smallest value (~ 1.1) near the midplane, implying 
that convection has been induced most strongly there. 


4.3 FU Ori Outbursts 


The FU Ori outbursts have characteristic amplitudes of ~ 6 
magnitudes, and the durations of the outburst and the qui¬ 
escence are evaluated as ~ 10 2 years and ~ 10 3 years, re¬ 


spectively (Hartmann & Kenyon 1996). Some authors have 
applied the DIMs to attribute the outburst cycle to the 


thermal-viscous limit cycle (Kawazoe & Mineshige 1993 


Bell & Lin 1994). In the DIMs, a is a free parameter and 


is chosen so that the observed outburst cycle is reproduced. 
Choosing the same a both on the upper and lower branches 
generally produces a smaller amplitude and a smaller duty 
cycle (= the ratio of the outburst duration to the quies¬ 
cence duration) than the observations. Therefore, the upper 
branch of a high a and the lower branch of a low a are usu¬ 
ally combined. (This is because a thermal equilibrium curve 
generally shifts to lower right as a decreases.) Specifically, 
Phot needs to be larger than a CO oi by a factor of 10 to re¬ 
produce the observed amplitude and duty cycle. Here, cthot 
and a C ooi denotes a on the upper and lower branches, re¬ 
spectively, which is a common notation in the DIMs. Also, 
to explain the absolute value of the outburst duration, it is 
required that c^hot ~ 10 -3 . 
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By contrast, we have uniquely identified from the first 
principles simulations that Qhot ~ 0.14 (near the low E 
end) and crcooi 0.03 as well as Emax/Emin 1.6. 
Comparing with, for example, Figure 3 in |Bell Sz Lin| 
(1994), tthot/<^cooi ~ 5 is smaller by a factor of two while 
Emax/E m i n ~ 1.6 is an order of magnitude smaller. One 
of the consequences of the smaller surface density contrast 
across the bistability would be a smaller amplitude of the 
outbursts. Actually, if we supposed a thermal-viscous limit 
cycle on our thermal equilibrium curve (neglecting the subtle 
middle branch) in Figure [l] the mass accretion rate would 
be ~10 _5 M Q yr _1 in the outburst and ~10 -7 M Q yr -1 in the 
quiescence, the former being an order of magnitude smaller 
than that suggested by observations. Furthermore, the dura¬ 
tion of the outburst would be too short with our ahot ~ 0.14; 
it is two orders of magnitude larger than that chosen in the 
DIMs to reproduce the outburst lasting ~ 10 2 years. Thus, 
we may conclude that the thermal-viscous limit cycle alone 
would not explain the FU Ori outbursts. Note, however, that 
the above discussion is true provided that the turbulence is 
driven by MRI both on the upper and lower branches. As 
we saw in Figure [5j temperatures in the lower branch so¬ 
lutions are typically ~ 10 3 K. Therefore, the lower branch 
solutions may need to be reconsidered taking account of the 
non-ideal MHD effects as well as the stellar irradiation that 
could reduce the effects. 

Recently, another model for FU Ori Outbursts has been 
proposed ( Armitage et al.pOOl Zhu et ah||2009 Martin & 


Lubow 2011): Since temperatures in protoplanetary discs 


(except the inner region) would be too low for good cou¬ 
pling between gas and magnetic field, a dead zone is ex¬ 
pected where MRI turbulence is suppressed and the accret¬ 
ing gas piles up. Then, eventually, the gravitational instabil¬ 
ity would drive turbulence, which dissipates to heat the disc 
while the disc is cooled by radiation ( Gammie||2001| Shi & 


Chiang|2014 ). Such thermal equilibria may compose another 


solution branch, the “gravito-turbulence” branch, below the 
classical S-shaped thermal equilibrium curve on the E-T e ff 
plane. Then, a limit cycle switching the gravito-turbulence 
branch and the MRI-turbulence branch would become possi¬ 
ble (Martin & Lu bow] 2011). This gravo-magneto limit cycle 
is shown to be reasonably consistent with the observations 
of FU Ori outbursts (Zhu et al. 2009). 

We comment on the rightmost two solutions on the up¬ 
per branch (ws0856 and ws0857). In these solutions, the 
shearing box approximation may not be valid since the 
pressure scale height h p of the disc patch is comparable 
or larger than the distance from the star Ro (See Table 
[TJ. Nevertheless, it might be worth to note that the ra¬ 
tio of vertically-integrated radiation pressure to vertically- 
integrated gas pressure exceeds unity (~ 1.7 in ws0856 and 
~ 2.9 in ws0857). They imply that the inner region of proto¬ 
planetary discs can be radiation-pressure dominated when 
M > 10 _3 M©yr _1 , which is an order of magnitude higher 
than that of the FU Ori outbursts though. Also, we note that 
a long-term build-up of heat is observed in ws0857 (while it 
is not in ws0856), which in fact is indicated by a relatively 


small value of 5 log Q /S log Q + in Figure 12 (see discussion 


4.2). 


Finally, we note that the a values obtained in shearing 
box simulations are generally those of (local) quasi-steady 
states. On the other hand, a values cannot be deduced from 


observations of steady accretion discs, and they are usu¬ 
ally evaluated utilizing transient phenomena of the viscous 
timescale t v is at radius R, using the following relation (e.g. 
Kotko & Lasota 2012), 


tv 


R 2 


Vv 


r 2 q 


act 


(32) 


2 7 Q) is the viscosity coefficient and c s is 
One needs to keep the above difference 


where z/ v is(~ cm; 
the sound speed 
in mind when one compares a values between shearing box 


simulations and observations (King et al. 


al. 2012). Actually, Sorathia et al. (2012) showed that the 


2007 Cannizzo et 


a values in the initial transient can be orders of magnitude 
larger than those in the quasi-steady state in their global 
(unstratified, isothermal) simulations. 


5 SUMMARY 

Using 3D radiation MHD simulations with realistic mean 
opacities and EOS, we explored thermodynamics in the in¬ 
ner part of protoplanetary discs where MRI turbulence is 
expected. Basic thermal properties are similar to those of 
dwarf nova discs in Paper I in spite of a two-orders-of- 
magnitude smaller Q. The thermal equilibrium curve con¬ 
sists of the upper, lower, and middle branches. The upper 
(lower) branch corresponds to a hot (cool) and optically very 
(moderately) thick disc patch, respectively, while the middle 
branch is characterized by convective energy transport near 
the midplane. Convection is also the major energy trans¬ 
port mechanism near the low E end of the upper branch. 
There, convective motion is fast in terms of Mach numbers 
reaching > 0.01, which enhances both MRI turbulence and 
cooling, raising a up to 0.14 that other wise is ~ 0.03. The 
enhancement of a due to convection near the low E end of 
the upper branch was also found in Paper I, and thus seems 
robust regardless of U of the disc patch. 

We examined causes of the S-shape of the thermal equi¬ 
librium curve based on the dependence of T m id on T e s and 
the time and horizontally-averaged vertical profiles plotted 
on the T-p plane. Then, we discussed thermal stability of 
the equilibrium solutions based on their trajectories on the 
logQ + -logQ _ plane. We also compared our results with 


the DIMs used to explain FU Ori outbursts (Kawazoe & 


Mineshige 1993 Bell & Lin 1994). Although the thermal 


equilibrium curve in our results also exhibits bistability, the 
surface density contrast across the bistability is an order 
of magnitude smaller, and the stress-to-pressure ratios in 
both upper and lower branches are two orders of magnitude 
greater, than those favored in the DIMs. It therefore ap¬ 
pears likely that FU Ori outbursts are not due solely to a 
thermal-viscous limit cycle resulting from accretion driven 
by local MRI turbulence. Instead, the limit cycle switch¬ 
ing the gravito-turbulence branch and the MRI-turbulence 
branch may be working there ( Armitage et al.||2001 Zhu et 


aL||2009 Martin &; Lubow~||2011 ). 

We note some caveats in our results. We ignored the 
stellar irradiation, although it is major heating source in 


4 Since c s can be expressed as a function of R based on the steady 
a model (Shakura &; Sunyaev 1973]), a can be evaluated once the 
viscous timescale t v is and the radius R are given. 
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protoplanet ary discs, to concentrate on the intrinsic ther¬ 
modynamics of accretion discs. Also, non-ideal MHD effects, 
which we didn’t include, might be important in the lower 
(cool) branch solutions. We are planning to take account of 
those as well as the net vertical magnetic flux in the future 
works. 
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Table 1. List of runs. 


run 


So 

Teffo 

ao 

ho 

E 

Teff 

a 

^total 

N x 

Ny 

N z 

L 'x 

L 'y 

L 'z 

L x /h p 

hp/ Ro 

^therm 

A 

ws0857 

Ue 

367337 

19952 

0.0292 

7.29e+ll 

365870 

26982 

0.0250 

334485 

32 

64 

384 

1.0 

4.0 

12.0 

1.03 

1.20 

12.13 

8.25 

ws0856 

U 

230461 

16982 

0.0294 

6.49e+ll 

228452 

18022 

0.0162 

547172 

32 

64 

384 

0.7 

2.8 

8.4 

1.07 

0.72 

16.82 

5.95 

ws0814 

U 

101199 

12302 

0.0295 

6.39e+ll 

119438 

13658 

0.0264 

312873 

32 

64 

256 

0.6 

2.4 

4.8 

1.20 

0.54 

10.29 

8.74 

ws0822 

U 

81609 

11220 

0.0295 

4.80e+ll 

77389 

12703 

0.0363 

222433 

32 

64 

256 

0.5 

2.0 

4.0 

0.95 

0.43 

7.79 

11.55 

ws0805 

U 

61523 

10000 

0.0293 

4.44e+ll 

60080 

9663 

0.0256 

217544 

32 

64 

256 

0.5 

2.0 

4.0 

1.12 

0.34 

11.60 

8.62 

ws0803 

U 

33289 

7943 

0.0294 
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The characters in the second column represents the outcome: U 
= upper branch, M = middle branch, L = lower branch, c = 
runaway cooling, and e = runaway heating. The units of surface 
densities (Eo and E), effective temperatures (T e ff 0 and T e ff), 
height (ho), and thermal time (t t h) are, respectively, gem -2 , K, 
cm, and orbit. L ' x , L ' y , and L' z are the box size (normalized by 
ho), and N x , N y , and N z are the number of cells, in the x , y, 
and in z directions, respectively. A is the period of 
time-averaging used in the diagnostics, normalized by ttherm- 
The pressure scale height is computed as 
hp = f (ptherm) dz/2 max((ptherm))• 
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Figure 1. Time-averaged effective temperature T e ff vs. surface 
density E of runs that reached a thermal equilibrium (large cir¬ 
cles). Large upward (downward) triangles indicate solutions that 
showed runaway heating (cooling) in the end, respectively. The 
initial condition of each run (T e ff 0 ,So) is shown as a small cir¬ 
cle, which is connected to (T e ff, E) by a gray line. Small upward 
(downward) triangles indicate the initial conditions of runs that 
didn’t reach a thermal equilibrium and showed runaway heat¬ 
ing (cooling), respectively. Colors represent the time-averaged 
Rosseland-mean optical thickness f t otal- Gray thick curves are 
thermal equilibria produced by a DIM without convection (see 
Appendix in Paper I), for a = 0.1, 0.03, and 0.01, from the top to 
the bottom. The name of each run is printed next to its (T e ffcb ^o) 
point with a small (but scalable) font for reference. 
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Figure 2. Time and horizontally-averaged vertical profiles of 
radiative heat flux F~ ad (red), advective heat flux F~ dv (blue), 
and cumulative heating rate t (gray), for selected solutions. 
In the inset in each panel, the position of the solution on the S— 
T e ff plane (Figure^) is indicated by a red circle. In each panel, the 
fluxes are normalized by the value shown on the right axis. The 
purple curve shows the hydrodynamic Brunt-Vaisala frequency 
squared divided by the angular velocity squared, N 2 /Q 2 . The 
plasma beta is larger than unity at the heights between the two 
vertical dotted lines. 
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Figure 3. Snapshots of various quantities on an x-z plane (y = 0) 
for a radiative solution (upper, t = 107.7 orbits in ws0805) and 
a convective solution (lower, t = 73.7 orbits in ws0837). From 
the left to the right, specific entropy, vertical advective heat flux 
ev z , density p, gas temperature T, and the total stress w xy are 
shown. The x and z axes are normalized by the time-averaged 
pressure scale height h p (see Table Q. Note that images here do 
not include the entire vertical extent of the box, but instead are 
limited to the midplane regions. Movies of snapshots of these two 
simulations are available online. 
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Figure 4. Time and horizontally-averaged vertical profiles of 
density (green), gas pressure (black), magnetic pressure (red), 
and radiation pressure (blue) for the selected solutions in Figure[2] 
The unit of pressures is ergcm -3 and the axis for density (gem - © 
is on the right. The vertical dotted lines denote the heights where 
the Rosseland-mean optical depth from the top/bottom boundary 
is unity. Other notations are the same as in Figure [2] 


Figure 5. Time and horizontally-averaged vertical profiles of gas 
temperature (K) (black) and ionization fraction (green) for the 
selected solutions in Figure [2] The axis for the ionization fraction 
is on the right. Other notations are the same as in Figure^] 
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Figure 6. Time-averaged quantities as a function of surface den¬ 
sity E: (a) a , (b) vertically-integrated stress Wa^(xO.l) (circles) 
and vertically-integrated thermal pressure thermal (triangles), 
each divided by E 4 / 3 , (c) advective fraction in the energy trans¬ 
port /adv, and (d) Mach number of the convective motion M a dv 
Colors represent the time-averaged effective temperature T e ff; the 
reddish and yellowish symbols correspond to the upper branch so¬ 
lutions while dark-greenish and dark-bluish symbols correspond 
to the middle and the lower branch solutions, respectively. 
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Figure 7. Time-averaged vertically-integrated Reynolds stress 
vs. Maxwell stress (upper) and vertically-integrated total stress 
W xy vs. f dz (lower). The dotted line indicates the equality. 
Other notations are the same as in Figure [6] 
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Figure 8. Convergence metrics (a) Q z , (b) Q y , (c) a ma g, and 
(d) /mag (see text for the definitions) as a function of surface 
density S. The gray horizontal line in panels (a) and (b) indicates 
a criterion for the numerical convergence while that in panels (c) 
and (d) indicates a criterion for well-developed MRI turbulence, 
suggested by |Hawley, Gu an, & Krolik (201l|. Other notations are 
the same as in Figure [6] 
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Figure 9. Time-averaged midplane temperature T m id = 
(T) (z = 0) vs. effective temperature T e ff- The dotted line in¬ 
dicates the critical slope, T m ^ oc T^ ff . Other notations are the 
same as in Figure [6] 



© 0000 RAS, MNRAS 000 , 000-000 











16 S. Hirose 



Gas Temperature (K) 


Figure 10. Time and horizontally-averaged vertical structure 
((T) ( 2 ), ( p) (z)) plotted on the T—p plane for all runs that reached 
a thermal equilibrium. The midplane and the photosphere cor¬ 
respond to the top and the bottom of each curve, respectively. 
Curves are colored by the advective fraction / a( j v in panel (a) 
while the background color indicates the Rosseland-mean opacity 
k; r (T, p) cm 2 g -1 in panel (b) and the first generalized adiabatic 
exponent Ti(T, p) in panel (c). 






Figure 11. Time trajectory of (Q + (£), Q~ (t)) during the time 
period A (see equation |10| >, for selected solutions. The portion 
corresponding to the last 15 orbits is indicated by color red for 
wa0831 and ws0850, which finally showed runaway cooling and 
heating, respectively. The displayed ranges of and Q~ is al¬ 
ways two orders of magnitude. Other notations are the same as 
in Figure [2] 
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Figure 12. The ratio of standard deviation of logarithmic cool¬ 
ing rate logQ _ (t) to that of logarithmic heating rate logQ + (t) 
(see Figure ED- The upward (downward) triangles indicate the 
solutions that showed runaway heating (cooling) in the end, re¬ 
spectively. Other notations are the same as in Figure [6] 
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APPENDIX A: RADIATION DAMPING AS A 
HEAT SOURCE 

Figure |A1| shows time and horizontally-averaged vertical 
profiles of cooling and heating rates (see equation [l7|), 


total U') = (<ldiss) - <(V • «)p - Vw : P) , (Al) 
q+m P (z) = - <(V • v)p - Vv : P) , (A2) 

9T * dWs (^ L )’ (A3) 

oZiv( z ) = (^(( e + - E )^))’ (A4) 


which are differential versions of those shown in Figure 
[2] Therefore, here we confirm again thermal equilibrium, 
where the total cooling rate and the total heating rate 
match. We note that most of the dissipation occurs in¬ 
side the Rosseland-mean photospheres in all cases, which 


justifies the first equality in equation (13). The compres- 
sional heating rate should be zero when vertically- 

integrated if the fluid motion is adiabatic. However, it 
can be positive when the fluid motion is damped by ra¬ 
diative diffusion, which is the case here. For example, 
in ws0800, f q+ompdz ~ 0.096 f g+ tal dz. That is, about 
10 % of the total heating comes from the virtual dis¬ 
sipation associated with compressional heating. We also 
computed the vertical component of f qtompdz to find 
f (p) (d (v z ) /dz)dz ~ —0.019 f q+ ota \dz. This means that 
the compressional heating comes mostly from horizon¬ 
tal compressions (= f qtompdz — f (p) (d (v z ) /dz)dz ~ 
0.115 f q^tgjdz), presumably due to MRI turbulence and 
the spiral acoustic waves. Below we evaluate the radiation 
damping rate q+ ad of an acoustic wave (that represent the 
fluid motion) to show q+ ad ~ q4 mp , which will conhrm that 
the non-zero (positive) compressional heating comes from 
the radiation damping. 

First, we start from the dispersion relation given in 


equation (1) in Blaes, Hirose, & Krolik 2007]) to get the 


exponential damping rate of the amplitude of a linear, plane 
acoustic wave: 


c j(e + 4:E)c 2 + (ick 2 /3k Fp)^Ec 2 ^ 2 

uj(e + 4 E) + (ick 2 /3 kfp)^E 


(A5) 
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^ 2 7 2 
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= c t k 2 
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1 - 
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(ck 2 /3 kfp)^E 
l j(e + AE) 
(ck 2 /3 kfp) 4:E 
c t k(e + 4 E) 


1 - 


(ick 2 /3 kfp)4:E^ ^ 2 


i_i 

2 

Ct 


u(e + E) J 


(A6) 

(A7) 

(A8) 


where uj and k are the frequency and the wave number of the 
acoustic wave, respectively, kf the flux-mean opacity, a the 
isothermal sound speed, and ct the total sound speed. From 


(A5) to (A6), we take the slow diffusion limit, where the 


first term dominates the second term both in the numerator 


and the denominator in the right hand side of (A5). Then, 
we take uj 2 s= c^k 2 (the undamped acoustic wave) to lowest 


order, and substitute it in the right hand side of ( |A7| ) to get 
to (A8), where uj is expressed as a function of k. Then, the 


damping rate of the amplitude T is obtained as the imagi¬ 


nary part of uj: 

1 (ck 2 /3K F p)^E 


r = - 


2 (e + 4 E) 


1-4 


2 ck 2 


3kfP E 
ck 2 




6 kfP 


(e»£) 
(e < E), 


(A9) 


(A10) 


where the factor 1/2 in the right hand side of (A9) comes 
from Taylor expansion of square root of ( |A8| . Hereafter, we 
assume that e <C E, and the total pressure is replaced with 
gas pressure. Then, as derived in equation (24) in Blaes et al 


(2011), the acoustic radiative damping rate can be evaluated 
as 


f/rad 


(fyn 


pc 2 

dpm&: 


fr 


= ^u-4- 


Ti 


2 ckp 


Ti J 3kr p/kE 


Sp n 


P 


(All) 

(A12) 

(A13) 


where Sp ma ^ is the pressure amplitude of the acoustic wave. 
From (A12) to (A13), (A10) was substituted with the ap- 

we can compute the time- 


proximation kf ~ Rr. 

Using simulation results, 
averaged version of (A 13) as 


9r + ad(2) = TFT 1 - TFT 


<Ti> 


2 ck ( p ) (e) ( Sp, 


(ri) J 3 (k r p) /k (E) V P 


(A14) 


where the wave number k and the pressure fluctuation 
dpmax/p may be evaluated from a snapshot of a simulation. 
Here, we take ws0800 as an example, and in that case, they 
are evaluated from Figure |A2| as 


k ^ 


27T 

(box width in x) 


Sp max 

P 


0.13. 


27T 

1.6 x 10 10 (cm) 


= 3.9 x 10 -lo (cm -1 
(A15) 


(A16) 


As shown in the panel of ws0800 in Figure m the resul¬ 
tant radiation damping rate q+ ad (z) shows a good agree¬ 
ment with the compressional heating rate q + 0mp (4 for the 
heights where qtomp( z ) is consistently positive. Therefore, 
the consistent positive compressional heating is regarded as 
the acoustic radiation damping. 
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Figure Al. Time and horizontally-averaged vertical profiles 
of radiative cooling rate q~ ad (red), advective cooling rate q~ dv 
(blue), the compressional heating rate (dotted) q £ Q mpi and the 
total heating rate ^ tal (gray), for the selected solutions in Fig¬ 
ure [2] normalized by the value shown on the right axis in each 
panel. The green curve in the pane l of ws0800 shows the acoustic 
radiation damping rate q+ ad (|A14l, which may not be applicable 
in high altitudes where the slow diffusion limit and the ignorance 
of magnetic field are invalid. Other notations are the same as in 
Figure [2] 
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Figure A2. Fluctuation of gas pressure relative to its horizontal 
average on the midplane at t = 202.0 orbits in the lower branch 
solution, ws0800. Axes are normalized by the time-averaged pres¬ 
sure scale height h p . 
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